module IO_mod

  use netcdf
  use namelist_mod
  use variables_mod

  implicit none

contains

subroutine read_nc_file()
  
  integer ncid, lat_dimid, lon_dimid
  integer start(2), count(2)
  integer lat_varid, lon_varid, u_varid, v_varid
  integer :: ndims_in, nvars_in, ngatts_in, unlimdimid_in
  integer :: status

  print*, 'Reading netCDF file: ', file_in

  call check(nf90_open(file_in, nf90_nowrite, ncid))
  call check(nf90_inquire(ncid, ndims_in, nvars_in, ngatts_in, unlimdimid_in))
  
  call check(nf90_inq_dimid(ncid, "lat", lat_dimid))
  call check(nf90_inquire_dimension(ncid, lat_dimid, len=nlats_in))

  call check(nf90_inq_dimid(ncid, "lon", lon_dimid))
  call check(nf90_inquire_dimension(ncid, lon_dimid, len=nlons_in))

  print*, 'nlons_in, nlats_in = ', nlons_in, nlats_in

  nlons = nlons_in - 2
  nlats = nlats_in - 2
  Nwave = min(nlons, nlats) - 1
  
  call allocate_vars(nlons_in, nlats_in, Nwave, nlons, nlats)

  call check(nf90_inq_varid(ncid, "lat", lat_varid))
  call check(nf90_inq_varid(ncid, "lon", lon_varid))
  call check(nf90_get_var  (ncid, lat_varid , lats))
  call check(nf90_get_var  (ncid, lon_varid , lons))

  dlon = (lons(2) - lons(1)) * deg2rad
  dlat = (lats(2) - lats(1)) * deg2rad

  call check(nf90_inq_varid(ncid, "u", u_varid))
  call check(nf90_inq_varid(ncid, "v", v_varid))

 
  count = (/nlons_in, nlats_in/)
  start = (/1, 1/)
  
  ! do record = 1, ntimes
  ! start(ndims_in) = record
  call check(nf90_get_var(ncid, u_varid, u, start=start, count=count))
  call check(nf90_get_var(ncid, v_varid, v, start=start, count=count))
  ! end do
  
  call check(nf90_close(ncid))

end subroutine read_nc_file

subroutine output(power1, power2, wavelength, wavenumber)

  real(r8), intent(in) ::  power1(Nwave), power2(Nwave), &
                           wavelength(Nwave), wavenumber(Nwave)

  integer :: ncid
  integer wavenumber_dimid, power1_varid, power2_varid, wavelength_varid, wavenumber_varid

  call check(nf90_create(file_out, nf90_clobber, ncid))
  call check(nf90_def_dim(ncid, "wavenumber", Nwave, wavenumber_dimid))

  call check(nf90_def_var(ncid, "wavenumber"    , nf90_double, wavenumber_dimid, wavenumber_varid))
  call check(nf90_put_att(ncid, wavenumber_varid, "units"      , "radians m-1"  ))

  call check(nf90_def_var(ncid, "power_var1"    , nf90_double, wavenumber_dimid, power1_varid))
  call check(nf90_put_att(ncid, power1_varid    , "units"      , "J kg-1"))
  call check(nf90_def_var(ncid, "power_var2"    , nf90_double, wavenumber_dimid, power2_varid))
  call check(nf90_put_att(ncid, power2_varid    , "units"      , "J kg-1"))
  call check(nf90_def_var(ncid, "wavelength"    , nf90_double, wavenumber_dimid, wavelength_varid))
  call check(nf90_put_att(ncid, wavelength_varid, "units"      , "km"    ))
  call check(nf90_enddef(ncid))
 
  call check(nf90_put_var(ncid, wavenumber_varid, wavenumber))
  call check(nf90_put_var(ncid, power1_varid  , power1))
  call check(nf90_put_var(ncid, power2_varid  , power2))
  call check(nf90_put_var(ncid, wavelength_varid, wavelength))
  call check(nf90_close(ncid))

end subroutine output

subroutine check(status)

  integer, intent(in):: status

  if (status /= nf90_noerr) then
    print*, trim(nf90_strerror(status))
    stop "stoped"
  end if

end subroutine check

end module IO_mod